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Abstract 

We present an algorithm for computing the determinant of an in- 
teger matrix A. The algorithm is introspective in the sense that it 
uses several distinct algorithms that run in a concurrent manner. Dur- 
ing the course of the algorithm partial results coming from distinct 
methods can be combined. Then, depending on the current running 
time of each method, the algorithm can emphasize a particular vari- 
ant. With the use of very fast modular routines for linear algebra, our 
implementation is an order of magnitude faster than other existing 
implementations. Moreover, we prove that the expected complexity 
of our algorithm is only 0(rt 3 log 2 ' 5 (ra||A||) ) bit operations in the case 
of random dense matrices, where n is the dimension and \\A\\ is the 
largest entry in the absolute value of the matrix. 

1 Introduction 

One has many alternatives to compute the determinant of an integer matrix. 
Over a field, the computation of the determinant is tied to that of matrix 
multiplication via block recursive matrix factorizations [19]. On the one 
hand, over the integers, a naive approach would induce a coefficient growth 
that would render the algorithm not even polynomial. On the other hand, 
over finite fields, one can nowadays reach the speed of numerical routines 
[12]. 

Therefore, the classical approach over the integers is to reduce the com- 
putation modulo some primes of constant size and to recover the integer 



determinant from the modular computations. For this, at least two variants 
are possible: Chinese remaindering and p-adic lifting. 

The first variant requires either a good a priori bound on the size of the 
determinant or an early termination probabilistic argument [13, §4.2]. It 
thus achieves an output dependant bit complexity of 0( log (| det (A) |) (n w + n 2 log (||-A||)) ) 
where to is the exponent of matrix multiplication 1 . Of course, with the co- 
efficient growth, the determinant size can be as large as O (nlog(n\\A\\)) 
(Hadamard's bound) thus giving a large worst case complexity. The algo- 
rithm is Monte Carlo type, its deterministic (always correct) version exists 
and has the complexity of O ((n w+1 + n 3 log(||A||)) log(n||A||)) bit opera- 
tions. 

The second variant uses system solving and p-adic lifting [6] to get a 
potentially large factor of the determinant with a O (n 3 log 2 (n|| J 4||)) bit 
complexity. Indeed, every integer matrix is unimodularly equivalent to a 
diagonal matrix S equal to diag (s±, . . . , s n ), where divides Sj+i. This 
means that there exist integer matrices U,V with det U, det V = ±1, such 
that A = USV. The Sj are called the invariant factors of A. In the presence 
of several matrices we will also use the notation Si(A). Solving a linear 
system with a random right hand side reveals s n as the common denominator 
of the solution vector entries with high probability, see [24, 1]. 

The idea of [1] is thus to combine both approaches, i.e. to approximate 
the determinant by system solving and recover only the remaining part 
(det(^4)/s n ) via Chinese remaindering. The Monte Carlo version of Chinese 
remaindering leads to an algorithm with the expected output-dependant 
bit complexity of C»(n w log (|^D|) + n 3 log 2 (n|| A\\)) . We use the notion 



of the expected complexity to emphasize that it requires E [log [i^jj to 
be O(l), where s n is the computed factor of s n and E denotes the expected 
value computed over all algorithm instances for a given matrix A. 

Then G. Villard remarked that at most O ^-y/log (| det (A) \ )j invariant 
factors can be distinct and that in some propitious cases we can expect 
that only the last O (log(n)) of those are nontrivial [17]. This remark, to- 
gether with a preconditioned p-adic solving to compute the i-th invariant 
factor lead to a O ^n 2+ t log 1 ' 5 (n|| J 4||) log°' 5 (n)^ worst case Monte Carlo 
algorithm. Without fast matrix multiplication, the complexity of the al- 
gorithm becomes O (n 3 ' 5 log 2 ' 5 (n|| A\\) log 2 (n)) . The expected number of 
invariant factors for a set of matrices with entries chosen randomly and uni- 



1 the value of u is 3 for the classical algorithm, and 2.375477 for the Coppersmith- 
Winograd method, see [4] 
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formly from the set of consecutive integers {0,1 A — 1} can be proven 

to be 0(log(n)). Thus, we can say that the expected complexity of the 
algorithm is O (n 3 log 2 (n||^4||) log(n) log^(n)) . Here, the term expected is 
used in a slightly different context than in algorithm [1] and describes the 
complexity in the case where the matrix has a propitious property i.e., the 
small number of invariant factors. 

In this paper we will prefer to use the notion of the expected rather than 
average complexity. Formally, to compute the average complexity we have 
to average the running time of the algorithm over all input and argorithm 
instances. The common approach is thus to compute the expected outputs 
of the subroutines and use them in the complexity analysis. This allows us to 
deal easily with complex algorithms with many calls to subroutines which 
depend on randomization. The two approaches are equivalent when the 
dependency on the expected value is linear, which is often the case. However, 
we can imagine more complex cases of adaptive algorithms where the relation 
between average and expected complexity is not obvious. Nevertheless, we 
believe that the evaluation of the expected complexity gives a meaningful 
description of the algorithm. We emphasize the fact, that the propitious 
input for which the analysis is valid can often be quickly detected at runtime. 

Note that the actual best worst case complexity algorithm for dense ma- 
trices is 0~ (n 2 - 7 log(||^4||)), which is 0~ (n 3 - 2 log(||^4||)) without fast ma- 
trix multiplication, by [21]. We use the notion 0~(N a log(||A||)), which 
is equivalent to 0(N a log^(N) log(||^4||)) with some (3 > 0. Unfortunately, 
these last two worst case complexity algorithms, though asymptotically bet- 
ter than [17], are not the fastest for the generic case or for the actually 
attainable matrix sizes. The best expected complexity algorithm is the 
Las Vegas algorithm of Storjohann [26] which uses an expected number of 
O (n w log (n 1 1 A\\) log 2 (n)) bit operations. In section 5 we compare the per- 
formance of this algorithm (for both certified and not certified variants) to 
ours, based on the experimental results of [27]. 

In this paper, we propose a new way to extend the idea of [25, 28] to get 
the last consecutive invariant factors with high probability in section 3.2. 
Then we combine this with the scheme of [1]. 

This combination is made in an adaptive way. This means that the algo- 
rithm will choose the adequate variant at run-time, depending on discovered 
properties of its input. More precisely, in section 4, we propose an algorithm 
which uses timings of its first part to choose the best termination. This par- 
ticular kind of adaptation was introduced in [23] as introspective; here we 
use the more specific definition of [5] . 
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In section 4.2 we prove that the expected complexity of our algorithm is 



0(n 3 log 2 (n||A||)0oiH) 

bit operations in the case of dense matrices, gaining a log L5 (n) factor com- 
pared to [17]. 

Moreover, we are able to detect the worst cases during the course of the 
algorithm and switch to the asymptotically fastest method. In general this 
last switch is not required and we show in section 5 that when used with 
the very fast modular routines of [9, 12] and the LinBox library [10], our 
algorithm can be an order of magnitude faster than other existing imple- 
mentations. 

A preliminary version of this paper was presented in the Transgressive 
Computing 2006 conference [14]. Here we give better asymptotic results 
for the dense case, adapt our algorithm to the sparse case and give more 
experimental evidences. 

2 Base Algorithms and Procedures 

In this section we present the procedures in more detail and describe their 
probabilistic behavior. We start by a brief description of the properties of 
the Chinese Remaindering loop (CRA) with early termination (ET) (see 
[7]), then proceed with the LargestlnvariantFactor algorithm to compute s n 
(see [1, 17, 25]). We end the section with a summary of ideas of Abbott et 
al. [1], Eberly et al. and Saunders et al. [25]. 

2.1 Output dependant Chinese Remaindering Loop (CRA) 

CRA is a procedure based on the Chinese remainder theorem. Determi- 
nants are computed modulo several primes Then the determinant is 
reconstructed modulo po- ■ -pt in the symmetric range via the Chinese re- 
construction. The integer value of the determinant is thus computed as soon 
as the product of pi exceeds 2| det (A) |. We know that the product is suffi- 
ciently big if it exceeds some upper bound on this value or, probabilistically, 
if the reconstructed value remains identical for several successive additions 
of modular determinants. The principle of this early termination (ET) is 
thus to stop the reconstruction before reaching the upper bound, as soon as 
the determinant remains the same for several steps [7]. 
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Algorithm 2.1 is an outline of a procedure to compute the determinant 
using CRA loops with early termination, correctly with probability 1 — e. 
We start with a lemma. 

Lemma 2.1. Let H be an upper bound for the determinant (e.g. H can 
be the Hadamard's bound: \det(A) \ < (y/n\\A\\ ) n ). Suppose that distinct 
primes pi greater than I > are randomly sampled from a set P with \P\ > 
2 [log; (H)~\ . Let r t be the value of the determinant modulo po - • -pt computed 
in the symmetric range. We have: 

'riog,(|det(A)|)l if det (A) ^0 
i/det(A)=0' 



(i) r t = del {A), if t > N = 



(ii) if rt 7^ det (A), then there are at most R = [log; primes 
Pt+i such that r t = det (A) mod po ■ ■ ■ PtPt+i; 

(Hi) if r t = r t +i = ■■■ = r t+k and ^p^j^^p^ < e, where R' = 
r io S/Sil' then P(r t ^ det (A)) <e. 

(iv) if r t = r t+1 = ■■■ = r t+k and k > \ lo^J^lg^H)^ ' where P ' = 
\P\ - [log; (H)] , then V (r t / det (A)) < e. 

Proof. For (i), notice that -[S^J < r t < [2^2*]. Then r t = det (A) 
as soon as po ■ ■ -p t > 2| det (A) \ . With / being the lower bound for pi this 
reduces to t > [log; | det (^4) |] when det (A) / 0. 

For (ii), we observe that det (A) = r t + Kpo . . .p t and it suffices to estimate 
the number of primes greater than I dividing K. 

For (iii) we notice that k primes dividing K are to be chosen with the prob- 

( R ) 

ability ,| P |_^/ +1) . . Applying the bound R' for R leads to the result. 

k 

For (iv) we notice that the latter is bounded by ^-p^ since R' < [log; (^-)] < 
\P'\. Solving for k the inequality (^jt^J < e gives the result. □ 

The two last points of the theorem give the stopping condition for early 
termination. The condition (iii) can be computed on-the-fly (as in Algorithm 
2.1). As a default value and for simplicity (iv) can also be used. 

To compute the modular determinant in algorithm 2.1 we use the LU 
factorization modulo pi. Its complexity is O (n w + n 2 log(||^4||)) . 

Early termination is particularly useful in the case when the computed 
determinant is much smaller than the a priori bound. The running time of 
this procedure is output dependant. 
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Algorithm 2.1 Early Terminated CRA 
Require: n x n integer matrix A. 
Require: < e < 1. 

Require: H - Hadamard's bound (H = (^/n\\A\\) n ) 

Require: / > 0, a set P of random primes greater than /, \P\ > 2[log ; (H)~\ . 
Ensure: The integer determinant of A, correct with probability at least 
1-e. 

1: i = 0; 

2: repeat 

3: Choose uniformly and randomly a prime pi from the set P; 

4: P = P\{ Pl } 

5: Compute det (A) mod pf, 

6: Reconstruct rj, the determinant modulo Po - ■ -pf, //by Chinese 
remaindering 

7: k = max{t : r^ t = ■■■ = n}; R' = [log, H+N 1 ; 

L 1 1 l} ' I toi p pi...pi-k 1 ' 

8: Increment i; 



2.2 Largest Invariant Factor 

A method to compute s n for integer matrices was first stated by V. Pan 
[24] and later in the form of the LargestlnvariantFactor procedure (LIF) 
in [1, 17, 7, 25]. The idea is to obtain a divisor of s n by computing a 
rational solution of the linear systems Ax = b. If b is chosen uniformly 
and randomly from a sufficiently large set of contiguous integers, then the 
computed divisor can be as close as possible to s n with high probability. 
Indeed, with A = USV, we can equivalently solve SVx = U~ 1 b for y = Vx, 
and then solve for x. As U and V are unimodular, the least common multiple 
of the denominators of x and y, d(x) and d(y) satisfies d(x) = d(y)\s n . 

Thus, solving Ax = b enables us to get s n with high probability. The cost 
of solving using Dixon p-adic lifting [6] is O (n 3 log 2 (n 1 1 A \ \ ) + n log 2 ( 1 1 b\ \ )) as 
stated by [22]. 

The algorithm takes as input parameters (3 and r which are used to 
control the probability of correctness; r is the number of successive solvings 
and (5 is the size of the set from which the values of a random vector b are 
chosen, i.e. a bound for ||6||. With each system solving, the output s n of the 
algorithm is updated as the 1cm of the current solution denominator d(x) 
and the result obtained so far. 
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The following theorem characterizes the probabilistic behavior of the LIF 
procedure. 

Theorem 2.2. Let A be a nx n matrix, H its Hadamard's bound, r and (5 
be defined as above. Then the output s n of Algorithm LargestlnvariantFactor 
of [1] is characterized by the following properties. 

(i) If r = 1, p is a prime, I > 1, then V (Vlf^) < \ \^\\ 

(ii) ifr = 2,P= [log (H)] then E (log (§*)) = O (1) ; 
(Hi) if r = 2, (3 = 6+ |~21og (H)~\ then s n = s n with probability at least 1/3; 
(iv) ifr= [2 log (log (#))!, (3>2 tten E (log (V)) = O (1) ; 



(v) if r = [log (log (H)) + log(i)], 2 | j3 and (3 > 2 then s n = s n with 
probability at least 1 — e; 

Proof. The proofs of (i) and (iv) are in [l][Thm. 2, Lem. 2]. The proof 
of (hi) is in [17][Thm. 2.1]. To prove (ii) we adapt the proof of (hi). The 
expected value of the under-approximation of s n is bounded by the formula 



E E nwfMi) . 

p\s n k=l \P P / 



where the sum is taken over all primes dividing s n . As [~-^] is bounded by 
i 

(3 



4 + \ this can be further expressed as 



oo oo Uog p (sn)J 

p prime k=l ' p\s„ k=l ' p|s„ fe=l 

E lo s (p) ^231 + a E lo s (p) + ~m E lo § (p) lo §p 

P prime p | s „ p | s „ 

1.78 + ^f») + l^)<5€0(l). 

To prove (v) we slightly modify the proof of (iv) in the following manner. 
From (i) we notice that for every prime p dividing s n , the probability that 
it divides the missed part of s n satisfies: 

Sn\ . I 1 



V p ^ < 
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As there are at most log (H) such primes, we get 

V (s n = s n ) > 1-log (H) (l/2) r > 1-log (H) 2-MiogW)-iog(|) = ^ {H) _L^. 

log (H) 

□ 

Remark 2.3. Theorem 2.2 enables us to produce a LIF procedure, which 
gives an output s close to s n with the time complexity O ^n 3 (log (n) + log (||^4||)) 2 ^ 
(see (ii)). 



2.3 Abbott-Bronstein-Mulders, Saunders- Wan and Eberly- 
Giesbrecht-Villard ideas 

Now, the idea of [1] is to combine both the Chinese remainder and the LIF 
approach. Indeed, one can first compute s n and then reconstruct only the 
remaining factors of the determinant by reconstructing det (A) / s n . The ex- 
pected complexity of this algorithm is O (n w log (| det(^4)/s n |) + n 3 log 2 (n||A||)) 
which is unfortunately 0~ (n w+1 ) in the worst case. 

Now Saunders and Wan [25, 28] proposed a way to compute not only s n 
but also s ra _i (which they call a bonus) in order to reduce the size of the 
remaining factors det(A)/ (s n s n _i). The complexity doesn't change. 

Then, Eberly, Giesbrecht and Villard have shown that for the dense case 
the expected number of non trivial invariant factors is small, namely less 
than 3 [log^ (n)] + 29 if the entries of the matrix are chosen uniformly and 
randomly in a set of A consecutive integers [17]. As they also give a way 
to compute any Sj, this leads to an algorithm with the expected complexity 
O (n 3 log 2 (n||,4||) log (n) log A (n)). 

Our analysis yields that the bound on the expected number of invariant 
factors for 

a random dense matrix can be refined as O (log ' 5 (n)). 

Then our idea is to extend the method of Saunders and Wan to get 
the last invariant factors of A slightly faster than by [17]. Moreover, we 
will show in the following sections that we are able to build an adaptive 
algorithm solving a minimal number of systems. 

The analysis also yields that it should be possible to change a log (n) 
factor in the expected complexity of [17] to a log(log(n)). This would 
require a small modification in the algorithm and a careful analysis. As- 
suming that the number of invariant factors is the expected i.e. it equals 
N = O (log (n)), we can verify the hypothesis by computing the (n— N— l)th 
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factor. If it is trivial, the binary search is done among O (log (n)) el- 
ements and there are only O (log (n)) factors to compute, which allows 
to lessen the probability of correctness of each OIF procedure. Thus, in 
the propitious case, the expected complexity of the algorithm would be 
O (n 3 log 2 (n||^4||) log A (n) log 2 (log(n))) . However, this cannot ce considered 
as the average complexity in the ordinary sense since we do not average over 
all possible inputs in the analysis. 



3 Computing the product of O (log(n)) last invari- 
ant factors 

3.1 On the number of invariant factors 

The result in [17] says that a n x n matrix with entries chosen randomly and 
uniformly from a set of size A has the expected number of invariant factors 
bounded by 3[log A (n)] + 29. In search for some sharpening of this result 
we prove the following theorems. 

Theorem 3.1. Let A be an nx n matrix with entries chosen randomly and 
uniformly from the set of contiguous integers {— [^J ... \^]}- Let p be a 
prime. The expected number of non-trivial invariant factors of A divisible 
by p is at most 4- 

Theorem 3.2. Let A be an nx n matrix with entries chosen randomly and 
uniformly from the set { — . . . The expected number of non trivial 

+ 3. 



invariant factors of A is at most 



V 2l °Sx (n) 



In order to prove the theorems stated above, we start with the following 
lemmas. 

( 1 A \ J 

Lemma 3.3. If j >1 the sum ( t[ — 1 j over primes p can be upper 

8<p<>S p ' 

bounded by 

Proof. We will consider separately the primes from the interval t^tt < P < 
t^, k = 0, 1, . . . k max . The value of k max is computed from the condition 
p > 8 and is equal to [log (A)] — 4. For the fcth interval [^] is less than or 
equal to 2 k+1 . In each interval there are at most [^zl numbers and 
at most P r i mes: if in the interval there are more than 3 odd numbers, 
at least one of them is divisible by 3 and is therefore composite. For this to 
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happen it is enough that \ 2 k m ax+5 ~\ — ^ which is the case. We may therefore 
calculate: 

n , sj riog(A)l-4 k+l , riog(A)l-4 /9 fc+ivj-l 1 riog(A)l-4 

E Qr^i < E 2^(V) - E V ) =«ft E(^ +1 ) 

8<p<A v 1 7 fc=0 v 7 fe=0 v 7 fc=0 

/flog(A)l-4 , NJ - 

□ 

Remark 3.4. For X = 2 l , k can be allowed from up to I — 3, instead of 
[log (A)] — 4 and we can include more primes in the sum. As a result we 

obtain an inequality Y,4< P <2' l)' ^ (i)'- 

Lemma 3.5. Xei Abeakxn, k<n integer matrix with entries chosen 
uniformly and randomly from the set {— |_^J . . . \^]} ■ The probability that 
rankp(A), the rank modulo p of A, is j, < j < k is less than or equal to 

j— 1 / -. \ maxfc— j— 1,0 

V (rank p (^) = j) < J[(l - a^) ■ p^k-il . j (1 + 0... (3 k ~ j ) 

w / lerea = _E_ L A±lj and/3= _i_ L A±lj. 

The proof of the lemma is given in the appendix A.l. 
Proof. (Theorem 3.1) 

The idea of the proof is similar to that of [17][Thm. 6.2]. 

For k > j let MDep fc (p, j) denote the event that the first k columns of 
A mod p have rank at most k — j over Z p . By Ij (p) we denote the event, 
that at least j invariant factors of A are divisible by p. This implies that the 
first columns of A have rank at most n — j mod p, or that MDep n _ J+fe (p, k) 
has occurred for all k = . . . j. This proves in particular that V (Ij (p)) < 
P (MDep n (p,j)). 

In order to compute the probability V (MDep s (p,j)), s > j we notice 
that it is less than or equal to 

s 

V (MDep, (p,j)) < V (MDep, (p,j))+ ^ V (MDe Pfc (p,j) A -< MDep fc _ 1 (p,j)) 

k=j+l 
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Surely, MDep^ (p, j) means that the first j columns of A are mod p, 
and consequently the probability is less than or equal to ftp 1 , where the value 
@P = ATI l"^p^"l ^ s a bound on the probability that an entry of the matrix 
is determined modulo p and is set to (A + l) -1 if p > A + 1 or less than or 
equal in the case p < A + 1. 

We are now going to find V (MDep fc (p,j) A -iMDep^ {p,j)) for k > j. 
Since the event MDep k _ 1 (p,j) did not occur, A k -i has rank modulo p at 
least (k — j) and of course at most {k — 1). For MDep^j to occur it must be 
exactly (k — j). This means that we can rewrite V (MDep fc (p,j) A -> MDep fc _ 1 
as 

V (MDe Pfe (p,j) | rank p (A fc _i) = k - j) ■ V (rank p (A k ^) = k - j) , 

where rank p (A k -i) denotes the rank modulo p of submatrix A k -\ of A, 
which consists of its first (k — 1) columns. 

Since the rank modulo p of A k _i is equal to k — j, there exists a set of 
k — j rows Lk-j which has full rank mod p. This means that we can choose 
k — j entries of the kth column randomly but the remaining n — k + j entries 
will be determined modulo p. This leads to an inequality 

P(MDep fc (p,j) | v a nk p (A k ^) = k-j)<(3;- k+ \ 

By Lemma 3.5 we have 7 3 (rank(Afc_ 1 = k—j) < (rrg - ) p( n - k +i)U 
Finally, we get 

V (MDe Pfc (p, j) A -i MDep k _ 1 (p,j)) < {j^j ^~ k+j)j (2) 

and 

(3) 

The expected number of invariant factor divisible by p < A verifies: 



E ( p & (p)) - p (p))) = E p & (p )) ^ E MDe p- (p- ■?') 



P-l/ VP + 1/ (p+l) J -2i 
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The latter is decreasing in p and therefore less than its value at p = 2, which 
is lower than 3.46. 

For p > A + 1 the result is even sharper: 

p ™ W ) - ^ (,))) , g ( i 4t)'" 1 Q)' x^t * i^tt^ 

the latter being lower than 1.18 for A > 1. □ 

Proof. (Theorem 3.2) 

In addition to MDep fc (p, j) introduced earlier, let Dep fc denote an event 
that the first k columns of A are linearly dependent (over rationals) and 
MDep fc (j), an event that either of MDep fc (p,j) occurred. 

Recall from [17, §6] that 

P(De Pl )<(A + l)- n 

^(DepfcA-Dep^) <V{Bep k | -Dep^) < (A + l)-^*" 1 . 

This gives V(Dep k ) < (A+ \ )n H h (A+1) i_ fc+1 which is less than (A+1) l_ fc+1 A j 1 • 

As in the previous proof, the probability that the number of non trivial 
invariant factors is at least j (event Ij) is lower than V (MDep n _j +k (k) V Dep n _ J+1 ) 
for all k = . . . j. The latter can be transformed to V ((MDep n _ J+fc (A;) A -> Dep n _ J+1 ) V Dep n _ J+1 ) , 
and both V (MDep n _ J+fc (/c) A ->Dep n _ J+1 ) andP (Dep n _ J+1 ) can be treated 
separately. 

To compute V (MDep n _j +k (k) A -i Dep n _ J+1 ) we will sum V (MDep n _ J+jfc (p, fc)) 
over all possible primes. Since Dep n _ J+1 does not hold, there exists a 
(n — j + 1) x (n — j + 1) non-zero minor, and we have to sum over the 
primes which divide it. We will treat separately primes p < A + 1 and 
p > A + 1. Once again we set f3 p = for p < A + 1 and (3 P = for 
p > A. 



By (3) we have 



X>(MDep„- J+t (P,*)> < (^'"V^ + (^'"ft^ 



p<A 

+ 



8<p<A 
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This transforms to 



kj nk / -[ \ kj 2^ 



X>(MDep n ^ +fc (p, k)) < 3*" 1 (§) 3 3^-, + 2*" 1 (0 













^ k-1 













2 k - 1 



37 3 k - 1 V 3 / V 4 / 4fe - 1 
6 fc / 1 r A + l.^' 



6 fc_l V A + 1 1 p 1 

8<p<A+l v 1 



Thanks to Lemma 3.3, the sum ^ 8<p<A+1 ^XTT r^p^"~l) ° can ' 3e bounded 

by . 

For primes p > A + 1 we should estimate the number of primes di- 
viding the (n — j + l)th minor. By the Hadamard's bound (notice that 
Dep n _j +1 does not hold), the minors are bounded in absolute value by 

( n-j+l 
(n — j + 1) (^2^)) ) • Therefore the number of primes p > A + 1 

dividing the minor is at most ^ (log A+1 {n) + 2) . Summarizing, 

P ((MDep B _ i+Jfc (*) A -n De Pn _ j+1 ) V Dep B _ i+1 ) < J ^— ^ + (^J 
3\ fc_1 f\\ k i 3 fc /4\ fe_1 /lA^ 4 fc /6A fe_1 6 k f l\ k ^ 



n /, , n ^/^ A + 1 \ fe_1 1 (A + l) fc A ,_ fn _,- +11 

We can now compute the expected number of non trivial invariant fac- 
tors. 

Let us fix h = max(2, \y/2 log A+1 (n) | ). We have that in particular, 
V(Ij) is less than V((MDep n _ j+h (h) A -> Dep n _ J+1 ) V Dep n _ J+1 ). We can 
check that h? > log A+1 (n) + log A+1 (log A+1 (n) + 2) . This gives also (A + 
l) h2 > n (log A+1 (n) + 2) and 

n „ , , v /A + lA^ 1 (A + l) 2/l 1 
1 > g (log A+1 (n) + 2) ^_ j (A + 1)h ( h+ i) - 

The expected number of non trivial invariant factors is bounded by: 

h n 

(h) V Dep n _ i+1 ) A Dep n _ J+1 ) 

j=l j=h+l 
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which in turn is bounded by 



f " {2\ hj 3 211 - 1 fl\ hj 2 2h ~ 1 fZ\ h - 1 {l\ h] 3 h 
~\3J \AJ 4*-l + V5/ e' 1 - 1 V2 ' 



1 



+ 5 (WO + *)(—) (A+ l)''i (A + 1)* - 1 + Hi^ + V > ' ' ') 

33/1-1 /2\M' l + 1 ) 2 3/i ^ 1 /3\ ft_1 3 2?i 

" k+ (3 h -2 h ) 2 {3) + (2 h -iy {2) + \2) (3 h -l) 2 [3) 

/4\ ft_1 4? h /i\ h ( h + 1 ) fQ\ h ~ 1 Q h 2 h f\\ h ( h+1 ) 
+ \3J (4 ft - l) 2 \A ) + \5J 6 h -12 fc -l (,2/ 

n h , , ^/A + lA^ 1 (A + l) 2/l 1 A + l\ 2 1 

+ 2( 1 °gA + iH + 2)^- r J ((A + l)^-l) 2 (A + 1)*(*+d H~J ATT 

, *t u " /, , n n x/A + lA ft_1 (A + l) 2?i 1 
< h + /(n, A) + 1. 

where / (n,Aj < (3 fc_ 2 /, )a ( 3 J +(2^-1)^ UJ +UJ (3^-1)^ UJ 

/4\^-l 4 2fl /l'vM'H- 1 ) 1 /'e^- 1 6 h 2 h (l\ h ( h + 1 ) 1 /'A+IA 2 _J_ ^ 9 QO 
"TV3^ (4 h -l) 2 U/ + l5i 6 h_i 2 h_il2/' V A / A+l ^ Z db 

soon as A > 2. On the other hand /(n, 1) = (2^1) I — ^ which leads to 
the final result. 

n 



3.2 Extended Bonus Ideas 



In his thesis [28], Z. Wan introduces the idea of computing the penulti- 
mate invariant factor (i.e. s n -i) of A while computing s n using two system 
solvings. The additional cost is comparatively small, therefore s n _i is re- 
ferred to as a bonus. Here, we extend this idea to the computation of the 
(n — k + l)th factor with k solvings in the following manner: 

1. The (matrix) solution of AX = B, where B is a n x k multiple right 
hand side can be written as s^N where s n approximates s n {A) and 
the factors of N give some divisors of the last k invariant factors of A: 
see lemma 3.6. 
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2. We are actually only interested in getting the product of these invari- 
ant factors which we compute as the gcd of the determinants of two 
perturbed k X k matrix R\N and R2N. 

3. Then we show that repeating this solving twice with two distinct right- 
hand sides B\ and B2 is in general sufficient to remove those extra 
factors and to get a very fine approximation of the actual product of 
the last k invariants: see lemma 3.10. 



3.2.1 The last k invariant factors 

Let X be a (matrix) rational solution of the equation AX = B, where 
B = [bi],i = 1, . . . , k, is a random n x k matrix. Then the coordinates of 
X have a common denominator s n and we let N = [ni\,i = 1, . . . , k, denote 
the matrix of numerators of X. Thus, X = s~ x N and gcd (iVy, s n ) = 1. 

Following Wan, we notice that s n (A) A~ l is an integer matrix, the Smith 
form of which is equal to 



diag ( 



(A) S n (A) Sn (A) 



\s n (A) ' S n -i (A)''"' Si (A) 

Therefore, we may compute s n -k+i (A) when knowing s& (s n (A) A~ Y ). The 
trick is that the computation of A^ 1 is not required: we can perturb A^ 1 
by right multiplying it by B. Then, su (s n (A) A~ l B) is a multiple of 
Sk (s n (A) A~ r ). Instead of s n (A) A~ 1 B we would prefer to use s n A~ 1 B 
which is already computed and equal to N. 
The relation between A and N is as follows. 

Lemma 3.6. Let X = s~ x N , gcd (s n ,N) = 1 be a solution to the equation 
AX = B, where B is n x k matrix. Let R be a random k x n matrix. Then 



gcd (Si (N) , s n ) 



'"n-i+l 



(A) and 



gcd ( Si (RN) , S n ) 



i+ i(A),i = l...,k. 



Proof. The Smith forms of s n (A) A 1 B and N are connected by the relation 
s -4^-Si (iV) = Si (s n (A) A~ l B), i = l,...,k. Moreover, Si (N) is a factor of 

Sl (RN). We notice that ^a^Ub) ec l uals 3^rJ> and thus gcd( 8i (W.O 
is an (integer) factor of s n -i+i [A). Moreover, the under-approximation is 
solely due to the choice of B and R. □ 

Remark 3.7. Taking gcd(sj (RN) ,s n ) is necessary as - may be a ra- 
tional number. 
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3.2.2 Removing the undesired factors 

In fact we are interested in computing the product ir^ = s n s n -\ ■ ■ ■ s n _k+i (A) 
of the k biggest invariant factors of A. Then, following the idea of [1], we 
would like to reduce the computation of the determinant to the computation 
of de ~^ , where tt^ is a factor of tt^ that we have obtained. We can compute 

7Tfc as §n/ gcd (/ife (RN) , §n) , where = S1S2 ■ ■ ■ Sk is the product of the k 
smallest invariant factors. 

We will need a following technical lemma. Its proof is given in the 
appendix, see A. 5. 

Lemma 3.8. Let V be an k x n matrix, such that the Smith form of V is 
trival. Let M be an nxk matrix with entries chosen randomly and uniformly 
from the set {a, a + 1 . . . a + S — 1} ; the probability that p l < S divides the 
determinant det(VM) is at most ^. 

In the following lemmas we show that by repeating the choice of matrix 
B and R twice, we will omit only a finite number of bits in We start 
with a remark, which is a modification of [28, Lem. 5.17]. We ramaind that 
the order modulo p (ord p ) of a value is the expotent of the highest power of 
p dividing it. 

Remark 3.9. For every n x n matrix M there exist a k x n, k < n, matrix 
V with trivial Smith form, such that for any nxk matrix B: if the order 
modulo p ordp ( Mfc ^fffi ) is greater than I then also ord p (det (VB)) is greater 



than /. 



Lemma 3.10. Let A be annxn integer matrix and Bi (resp.Ri), i = 1,2 be 
nxk (resp.k x n), matrices with the entries uniformly and randomly chosen 
from the set S of S contiguous integers, k > 2. Denote by [i^ the product 
s\ . . . Sk of the k smallest invariant factors and by ir^ the product of the k 
biggest factors of A. Then for M = s n (A) A^ 1 

log ( gcd ( W (fliMB!), ti k (R 2 MB 2 ), s n (A) k ) J J G O (l)+0 fl^iE. 



Sn(A) k V J J J \ S 

where H is the Hadamard bound for A. 
Proof. First, notice that nk ^} h = — hrv. Therefore 
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the expected value is less than or equal 
L L l og(p)lV ord p — 



P\Sn(A) 



= E E log(p)P (ord p ( 
i p|s„(A) V V 



MM) 

gcd (MRiMB 1 ),MR2MB 2 ), s n (A) k ) 



MM) 



> I 



P|Sn(A) 



° rd P I Mfc (M) J ^ fc A 



i p| S „(A) \k=0 \ OTd P ^ ft(MB.) ) Z (I ~ K) 

Thanks to remark 3.9 we can link this probability to the probability that 
p l divides the determinant of VBi or R4U, for matrices V, U which have a 
trivial Smith form. 

We only consider p\s n (A). 

For p l < S Lemma 3.8 gives 

\^vf a ( MMBj) \ ( M R jMBj) \ \ 

r dp {-mmtj - k ordp vmmb-) * {l - k) ) * 
1 3 

{Bi : ord p (det (VBi)) > k) V {R { : ord p (det (RiU)) >k)<(l + l) — . 
fc=0 p 

Now the expected size of the under-estimation is less than or equal to 
log (2) (3 + p ((/ + l) 2 |) ^ + log (3) ^2 + p ((/ + 1) |) ^ 

^(<-S(l)>'-«(g(l)> 5 |: ff S'-«© : 

, 4 , 6 + 2 ,4 + ,U + 0, 7+ E,o gW ^±^ 

/•OO 

J 10 (x-l) 3 (x + lf 

which is O (1). 



coo — 97r 2 -I- 3fi-r 4 -I- Q 

< 8.51 + / log (x) T + dx < 8.51 + 11.97 
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For p l > S the probability V (p l \ det (M)) is less than V (pL lo § P ( 5 )J | det (M) 

and consequently can be bounded by 3min (j^, which is less than 
The expected size of the underestimation is 

fclog p (ff) 

,-/ . ,^1,,,,,,,/ - V^i 'if^M„„. ,/,, . ^M„„ 

P 



E E (^D 2 io g (P)(^=) < E^| M (^ 1 ^w + ^( fc ^ 

p| Sn (A)«=riog p (S)l VVO/ p|a„(A) 

+ 1 (k log p (H)) 3 ) < k log 2 (IT) | + l k log + log 2 (//)) < ^ 



This is O ^ fc lo s ( H ) ^ ; which gives the result. □ 

Another method to compute the product of some first invariant fac- 
tors of a rectangular matrix N would be to compute several minors of the 
matrix and to take the gcd of them. In our scheme we can therefore get rid of 
matrix R which would enable us to use a smaller bound on S = O (k log (H)) 
and still preserve a small error of estimation due to the choice of B. How- 
ever, it is difficult to judge the impact of choosing only a few minors (instead 
of all). An experimental evaluation whether for random A and random B 
the minors of N are sufficiently "randomly" distributed remains to be done. 



4 Introspective Algorithm 

Now we should incorporate Algorithm 2.1 and the ideas presented in sections 
2.2 and 3.2 in the form of an introspective algorithm. 

Indeed, we give a recipe for an auto-adaptive program that implements 
several algorithms of diverse space and time complexities for solving a par- 
ticular problem. The best path is chosen at run time, from a self-evaluation 
of the dynamic behavior (here we use timings) while processing a given in- 
stance of the problem. This kind of auto-adaptation is called introspective 
in [5]. In the following, CRA loop refers to Algorithm 2.1, slightly modi- 
fied to compute det (A) /K. If we re-run the CRA loop, we use the already 
computed modular determinants first whenever possible. 

Informally, the general idea of the introspective scheme is: 

1. Initialize the already computed factor K of the determinant to 1; 

2. Run fast FFLAS LU routines in the background to get several modular 
determinants di = det (A) mod pi. 
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3. From time to time try to early terminate the Chinese remainder re- 
construction of det (A) /K. 

4. In parallel or in sequential, solve random systems to get the last in- 
variant factors one after the other. 

5. Update K with these factors and loop back to step (2) until an early 
termination occurs or until the overall timing shows that the expected 
complexity is exceeded. 

6. In the latter exceptional case, switch to a better worst case complexity 
algorithm. 

More precisely, the full algorithm in shown on page 20. 

4.1 Introspectiveness: dynamic choice of the thresholds 

The introspective behavior of algorithm 4.1 depends paramountly on the 
number of system solvings and on the size of the random entries. 

The parameter i ma x controls the maximal total number of system solv- 
ings authorized before switching to a best worst-case complexity algorithm. 
The choice of i max has to be discussed in terms of the expected number of 
invariant factors of A. 

First, depending on the size of the set from which we are sampling the 
random right-hand sides, a minimum number of solvings is required to get 
a good probability of correctness. We thus define this to be i m in- 

In the dense case, the (ii) part of theorem 2.2 states that i m i n = 2 
is sufficient. Part (iv) part of theorem 2.2 prompts us to take i m %n = 
[2 log (log (H))~\ if we want to use a smaller (3. 

Then this number i m i n has also to be augmented if the expected number 
of non trivial invariant factors is higher. We thus set 

B(# factors (A))). 

In the dense case E (# factors (A)) is less than \y/2 log A (n)] + 3 as shown 
in theorem 3.2 . 

Now, random vectors are randomly sampled a set of size S. For a dense 
matrix A we need S = 13E (# factors (A)) 3 ([log (H)~\) 4 to get a good prob- 
ability of success as shown in theorem 2.2(h) and lemma 3.10. 

Additionally, (see lemma 3.10) we should ensure that irk is computed 
twice using different matrices B. We therefore introduce the variables kd one 
and k app which store respectively the number of factors computed at least 
twice (up to 0(1)) or once (thus only approximated). 
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Algorithm 4.1 Extended Bonus Determinant Algorithm 
Require: An integer n x n matrix A. 

Require: H - bound for det(-A) (can be the Hadamard's bound) 
Require: < e < 1, an error tolerance, S = 

13E (# factors (A)) 3 ( [log (H)] f , I > 1. 
Require: A stream S of numbers randomly chosen from the set of S con- 
tiguous integers. 

Require: A set P of random primes greater than /, \P\ > \2logi(H)~\ , P' = 
\P\-\o gl {H) 

Ensure: The integer determinant of A, correct with probability at least 
1 -e. 

1: k = log(l/e) / [log (i^y)l; see Lem. 2.1(iv) 

2: for i = 1 to k do 

3: run the CRA loop for det (A) ; //see Alg. 2.1 

4: if early terminated then Return determinant end if 
5: end for 

6: ijnax = ^raax (A) , imin = ^min / /see §4.1 

7: n = l;K= 1; 

8 : ^done = Oj ^app = 0; J ' = 0; 

9: while /c rfone < i maa; do 

10: i = kdone + lj 

11: while i < i m „ ? do 

12: Generate 6 t a random vector of dimension n from the stream S 1 ; 

13: Compute s n by solving Axf ^ = b^; //see Section 2.2 

14: if i = 1 then ;7fi = s n ; 

15: else 

16: N := s n X, where A = [xp^ =lv ..j; //see Section 3.2; 

17: Generate a random ixn matrix i2. 

19: end if 

20: A = lcm (7Tj, A); 7Tj = A; 

21: Resume CRA looping on d = det (A) /A for at most the time of 
one system solving; 

if early terminated then Return d ■ A; end if 
if i > i min then 
if 7Tj = 7fj_i then 
if i > k app then 

kdone = k a pp', k a pp = %\ j ' = j ' + 1 mod 2; break; 
else 

Resume CRA looping on d = det (A) / A for at most the time 
of (i max ~ i) system savings; 
if early terminated then Return d ■ A; 
else i = i max ; end if 
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end if 
end if 
end if 

i=i+l; 
end while 



4.2 Correctness and complexity 

Theorem 4.1. Algorithm 4-1 correctly computes the determinant with prob- 
ability 1 — e. 

Proof. Termination is possible only by the early terminated CRA loop or 
by the determinant algorithm used in the last step. The choice of k from 
theorem 2.1(iv) and the choice of the determinant algorithm from [20, 27] 
ensures that 1 — e probability is obtained. □ 

The following theorem gives the complexity of the algorithm. 

Theorem 4.2. The expected complexity of Algorithm 4-1 in the case of a 
dense matrix is 

O (> log (1/e) + n 3 (logn + log (p||)) 2 log ' 5 (n)) . 

The worst case complexity depends on the algorithm used in the last step. 

Proof. To analyze the complexity of the algorithm we will consider the com- 
plexity of each step. 

For a dense matrix A, with k denned as in the line 1, the complexity of 
initial CRA iterations is O (n w log (1/e)). The while loop is constructed in 
this way that we perform at most 2i max (see subsection 4.1 for the bound 
on imax) iterations, where log(||-B||) = O (log(n) log(log(||^4||))). Therefore 
the cost is O ^n 3 (log(n) + log ( 1 1 A||)) 2 y / log(n)^ . Considering the time limit, 
this is also the time of all CRA loop iterations. To compute tti we 

need ni^^ 2 bit operations. Then, the computation of the i x i de- 
terminant of RN by a deterministic algorithm (i.e, deterministic CRA) 
costs 0(i tJ (log(i) + n(log(n||i?|| • \\A\\ ■ ||-B||)))) bit operations, which for 
i = 2,..., imax with imax being O (log (n)) is 0~ (n) and thus negligible. 

With the expected number of invariant factors bounded by i ma x (see 
Thm.3.2), it is expected that the algorithm will return the result before the 
end of the while loop, provided that the under-estimation of TTi max is not 
too big. But by updating s n O (log ' 5 (n)) times and updating the product 
TTi max twice, it is expected that the overall under-estimation will be 0(1) (see 
Theorem 2.2 and Lemma 3.10), thus it is possible to recover it by several 
CRA loop iterations. □ 

For the last step for a dense matrix we propose the 0~ (n 3 2 log (||A||)) 
algorithm of Kaltofen [21] or 0~ (n u log ( 1 1 A\\)) algorithm of Storjohann [27]. 
We refer to [20] for a survey on complexity of determinant algorithms. 
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5 Experiments and Further Adaptivity 



5.1 Experimental results 

The described algorithm is implemented in the LinBox exact linear algebra 
library [10]. In a preliminary version i max is set to 2 or 1 and the switch 
in the last step is not implemented. This is however enough to evaluate the 
performance of the algorithm and to introduce further adaptive innovations. 

All experiments were performed on 1.3 GHz Intel Itanium2 processor 
with 128 GB (196 GB since September 2006) of memory disponsible. 

For a generic case of random dense matrices our observation is that 
the bound for the number of invariant factors is quite crude. Therefore 
the algorithm 4.1 is constructed in the way that minimizes the number of 
system solving to at most twice the actual number of invariant factors for a 
given matrix. Under the assumption that the approximations s n and ff j are 
sufficient, this leads to a quick solution. 

Indeed for random dense matrices, the algorithm nearly always stopped 
with early termination after one system solving. This together with fast 
underlying arithmetics of FFLAS [9] accounted for the superiority of our 
algorithm as seen in figure 1 and 2 where comparison of timings for different 
algorithms is presented. Notice, that our algorithm beats the uncertified (i.e. 
Monte Carlo type) version of the algorithm of [26] which claims currently 
the best theoretical complexity. This proves that adaptive approach is a 
powerful tool which allow us to construct the algorithms very fast in practice 

Thank to the introspective approach our algorithm can detect the cases 
when the number of invariant factors is small and equal to k < i max ■ One 
can therefore argue the complexity of our algorithm is in fact O ^n 3 (log (n) + log (||A||)) 2 k 
where k is the number of invariant factors. To test the performance of our 
algorithm to detect propitious cases we have run it on various sets of struc- 
tured and engineered matrices. The adaptive approach allowed us to obtain 
very good timings which motivates us to encourage the use of this algorithms 
in the situations which go further beyond the dense case. 

Figure 3 we present the results of the determinant computation for sparse 
matrices of N. Trefethen 2 . 

The results encouraged us to construct a sparse variant of our algorithm, 
which we shortly describe in Section 5.2. Figure 3 gives a comparison of 
the performance of sparse and dense variants. We used the sparse solver of 
[18]. Using the algorithm with the dense solver outperforms using the sparse 

2 http:/ /ljk. imag.fr/membres/Jean-Guillaume.Dumas/Matrices/Trefethen/ 
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NTL (Monte Carlo) ■ 

LU-CRT ■ 

[Strojohann-Giorgi-Olesh (Certified)] ■ 
[Strojohann-Giorgi-Olesh (Uncertified)] 
'""ima (Monte Carlo) 

Hybrid algorithm I 



A )96 



Figure 1: Comparison of our algorithm with other existing implementation. 
Tested on random dense matrices of the order 400 to 10000, with entries 
{-8,-7,. . . ,7,8} Using fast modular routines puts our algorithm several times 
ahead of the others. Scaling is logarithmic. 



NTL (Monte Carlo) 
Abbott's LIF 
Hybrid algorithm 




Figure 2: Comparison of our algorithm with early terminated Chinese re- 
maindering algorithm (LU) and the algorithm of Abbott et al. [1] (LIF). 
Tested on random dense matrices of the order 40 to 1000, with entries {- 
100,-99,. . . ,99,100}. When matrix size exceeds 80 the adaptive algorithm 
wins. Scaling is logarithmic. 
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Figure 3: Comparison of sparse and dense variants of our determinant algo- 
rithm for Trefethen's matrices. Scaling is logarithmic. 

solver by a factor of 3.3 to 2.3, and decreasing with the matrix size n. Thanks 
to the space-efficiency of the sparse algorithm we are able to compute the 
determinant for 20000 x 20000 matrix for which the dense solver requires to 
much memory. 

In figure 4 we compare the performance of dense and sparse variants of 
the algorithm with the CRA algorithm (sparse variant) for random sparse 
matrices. The matrices are very sparse (20 non-zero entries per row). To 
ensure that the determinant is non-zero we put 1 on the diagonal. Both dense 
and sparse variants of the algorithm have better running times than the 
CRA, which proves that we can detect propitious cases for sparse matrices. 
Furthermore, sparse variant is best for bigger matrices and again lets us 
solve the problem when the dense variant fails due to unsufficient memory. 

In Table 1 W6 give the timings for the algorithm with imax 

= 1 and 2. 

The algorithms was run on a set of specially engineered matrices which have 
the same Smith form as diag{l, 2 . . . n} and the number of invariant factors 
of about S. We notice that the algorithm with i m ax = 1 (which is in fact 
a slightly modified version of Abbott's algorithm [1]) runs better for small 
n. This motivated us to develop an even more adaptive approach, which we 
describe in Section 5.3. 

5.2 Sparse matrix case 

When trying to adapt our determinant algorithm to the sparse case, the 
immediate problem is the bound for the expected number of invariant fac- 
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Figure 4: Comparison of sparse and dense variants of our determinant al- 
gorithm with the CRA algorithm for random sparse matrices. Scaling is 
logarithmic. The running time of the CRA algorithm has been approxi- 
mated based on the timings for one iteration 



n 


i"max — 1 


^max — 2 


n 


i"max — 1 


i"max — 2 


100 


0.17 


0.22 


300 


5.65 


5.53 


120 


0.29 


0.33 


350 


9.76 


9.64 


140 


0.48 


0.55 


400 


14.99 


14.50 


160 


0.73 


0.78 


600 


57.21 


54.96 


180 


1.07 


1.16 


800 


154.74 


147.53 


200 


1.49 


1.51 


1000 


328.93 


309.61 


250 


2.92 


3.00 


2000 


3711.26 


3442.29 



Table 1: Comparison of the performance of Algorithm 4.1 with z max set to 
1 and 2 on engineered matrices. 
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tors. On can easily notice, that for a matrix with k non-zero entries per 
row, chosen uniformly and randomly from the set of S contiguous integers, 
the expected number of invariant factors divisible by p l can be bounded 
from below by [^rj and thus is linear with n. Thus, we cannot use the 
same argument to estimate the expected number of system solvings as in 
the dense case. 

One solution would be to consider the number of "big" invariant fac- 
tors, i.e. the number of invariant factors which are bigger than a certain 
parameter C. The parameter has to be chosen in a way, that the product 
of all smaller factors can be computed by modular CRA loop quicker than 
another rational solution of a system of equation. We could exploit here 
the difference in complexity between system solving and one modular rou- 
tine which is 0~ (n 3 ) to 0(n w ) (or 0(n L5 f2) to 0(Qn) in the case of sparse 
procedures). This could enable us to recover O (n 3 ~ w ) (0(n°' 5 )) bits of the 
determinant by running modular routines without exceeding the cost of one 
linear system solving. This adaptive solution is already implemented in the 
dense version of the algorithm, which motivated us to run it on potentially 
unsuitable matrices. When comparing the times for the CRA algorithm and 
our algorithm applied to sparse matrices (however, without exploiting their 
sparsity) we decided that there is a need for a "sparse" version of our algo- 
rithm, which will take into account the sparse structure of the matrices in 
the subroutines. 

In what follows we will shortly present the sparse counterparts of the 
subroutines used, give their complexities and discuss some modification of 
the parameters if needed. We assume that the cost of one matrix- vector 
product is Q, = O (n) . 

Instead of the dense LU, sparse elimination can be used in practice e.g. 
for extremely sparse matrices [11]. In general, black box method are pre- 
ferred. The idea is to precondition the matrix so that its characteristic 
polynomial equals its minimal polynomial [16, 2]; and then to compute the 
minimal polynomial via Wiedemann's algorithm [29]. The complexity of 
the sparse modular determinant computation is then O (Qn) [11, Table 4]. 
Adaptive solutions exist [15]. 

For solving a sparse system of linear equations the solver of [18] can be 
used. By similar reasoning as in [22], the cost of solving Ax = b for a sparse 
matrix A is that of O (n L5 log(n|| J 4||) + n 0,5 log(||6||)) matrix-vector prod- 
ucts and O (n 2 log(n||A||) (n - 5 + log(||A||)) + n 2 log(||6||) log(n||A||) + n log 2 (||6||)) 
additional arithmetic operations. 

If ||6|| is O (n 0,5 ) in size and f2 = O(n), this means that the complexity 
of computing s n is 0(n 2 log(n||A||)(n a5 + log(|| j4|| ))) bit operations. 
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Currently known sparse determinant algorithms that can be used in the 
worst-case step include the CRA loop (with the complexity O (f2nlog(| det(A)|))) 
and the algorithm of [17]. By moving to the sparse solver in [17] we can ob- 
tain an algorithm with the worst time complexity of 0(n 3 log 1 ' 5 (n|| J 4||) log (|| A\\) log 2 (n) ) . 

All in all, by moving to the sparse procedures, we obtain the algorithm 
with the complexity O (min (k (On 15 log(n 1 1 A \ \ ) + n 2 - 5 log(n \\A \ \ ) log ( \\A \ | )) , n 2 log(ra 1 1 A \ \ ) 
where k is the number of invariant factors. In the propitious case where k is 
smaller than O (\/n) we obtain an algorithm with the running time better 
than the currently known algorithms. 

5.3 More adaptivity 

We start with a simple remark. For every matrix, with each step, the size 
of s n -i decreases whilst the cost of its computation increases. In Table 1, 
this accounts for better performance of Abbott's algorithm, which computes 
only s n , in the case of small n. For bigger n calculating s n _i starts to pay 
out. The same pattern repeats in further iterations. 

The switch between winners in Table 1 can be explained by the fact that, 
in some situations, obtaining s n _j by Ll/-factorization (which costs 
the time of LU) outperforms system solving. Then, this also holds for all 
consecutive factors and the algorithm based on CRA wins. The condition 
can be checked a posteriori by approximating the time of LUs needed to 
compute the actual factor. We can therefore construct a condition that 
would allow us to turn to the CRA loop in the appropriate moment. This 
can be done by changing the condition in line 27 (tt^ = 7Tj_i) to 

/ Tij \ time ( solving) , 

if the primes used in the CRA loop are greater than I. This would result 
with a performance close to the best and yet flexible. 

If, to some extend, s n -i could be approximated a priori, this condition 
could be checked before its calculation. This would require a partial factor- 
ization of s n -i+i and probability considerations as in section 3.1 and [17]. 

6 Conclusions 

In this paper we have presented an algorithm computing the determinant 
of an integer matrix. In the dense case we proved that the expected com- 
plexity of our algorithm is 0(n 3 log 2 (n||^4||) log ' 5 (n) ) and depends mainly 
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on the cost of the system solving procedure used and the expected number 
of invariant factors. Our algorithm uses an introspective approach so that 
its actual expected complexity is only O ^n 3 (log(n) + log(||A||)) 2 k^j if the 
number k of invariant factors is smaller than a priori expected but greater 
than i m i n ; The actual running time can be even smaller, assuming that any 
under-estimation resulting from probabilistically correct procedures can be 
compensated sooner than expected. Moreover, the adaptive approach allows 
us to switch to the algorithm with best worst case complexity if it happens 
that the number of nontrivial invariant factors is unexpectedly large. This 
adaptivity, together with very fast modular routines, allows us to produce 
an algorithm, to our knowledge, faster by at least an order of magnitude 
than other implementations. 

Ways to further improve the running time are to reduce the number 
of iterations in the solvings or to group them in order to get some block 
iterations as is done e.g. in [3]. A modification to be tested, is to try to 
reconstruct s n with only some entries of the solution vector x = n/d. 

Parallelization can also be considered to further modify the algorithm. 
Of course, all the LU iterations in one CRA step can be done in parallel. 
An equivalently efficient way is to perform several p-adic liftings in parallel, 
but with less iterations [8]. There the issue is to perform an optimally 
distributed early termination. 
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A Properties of matrices with almost uniformly 
distributed entries 



In this appendix we present some probabilistic properties of matrices with 
entries almost uniformly distributed modulo p l , I £ Z. We consider the 
case, when the entries are randomly and uniformly chosen from a set of S 
contiguous integers S = {a, a + 1 . . . a + S — 1}, for any a. As the result, 
the probability that an entry is equal to a given d modulo p l is bounded as 
follows 

2_ L^j < v{x:x = d modp *)<J_r^. (4) 
bp 1 bp 1 

We set 

This special case of non-uniformly distributed random variables was widely 
considered in the thesis of Z. Wan (see [28]) for I = 1. In the following we 
will first consider the rank modulo p of a matrix under certain conditions 
(lemma A.l). Then we give the analogues of the theorems 5.9-5.15 of [28] 
in the case I > 1 (lemmas A. 2, A. 3, A. 4). This allows us to prove Theorem 
3.2 on the expected number of invariant factors and Theorem 3.10, which 
gives the expected size of over-approximation of m in the case of perturbed 
matrices. 

Lemma A.l. Let Abeakxn, k<n integer matrix with entries chosen 
uniformly and randomly form S. The probability that rank p (^4) ; the rank 
modulo p of A, is j, < j < k is less than or equal to 

j- 1 / -, \ max(fc-j'-l,0) 

V (rank p (^) = j) < JJ(1 - a^) ■ /^X^) • — _ {1 + 0... P^) 

i=o ^ ^ 

<p{n-m-j)(_L_j \ (6) 

where a = [f J an d @ = ^ [f J • 

Proof. The proof is inductive on k — j and j. For j = and k leqn the fact 
that rankp(A) = means that all the entries of A are zero modulo p, that is 

V(rank p (A) = 0) < (3 n \ 
the latter being less than [5 nk (^j+^j ■ 
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Now, denote by Ai the submatrix of A consisting of i first columns. For 
k = j we have 

V(rank p (Ak) = k) = V(r&nk p (Ak) = k | rank p (A£;_i = k — 1)) ■ V(rank p (Ak~i = k — 1) 
k 

= P(rank p (yli) = 1) P(rank p (^4j) = i | rank p (,4j_i = i - 1)). 

To compute / 9 (rank p ( J 4j) = i | rank p ( J 4j_i) = z — 1) we notice the fact 
that rank p (^4j_i) = i — 1 means that we can choose an (i — 1) x (i — 1) 
non-zero minor of A^±. This means that we can leave the choice of the 
corresponding i — 1 entries of the ith column free and only have to ensure 
that the remaining subvector of size n — i + 1 is not equal to some given 
vector. This gives 

P(rank p (^i) = i | rank p (^_i = i - 1) < (1 - a n ~ i+1 ) 

and in consequence 

k 

V(vank p (A k ) =k)< - a n ~ i+1 ). 

i=l 

Now, assume that for all (j, k) such that k — j < M the bound (6) holds. 
We consider V(r&nk p (AK) = J), where K — J = M > 0. We can rewrite: 

"P(rank p (^4ft-) = J) = "P(rank p (^4ft-) = J | rank p ( J 4^„i) = J) • "P(rank p ( J 4^„i) = J) 
+ V(rank p (A K ) = J | rank p (A^_i) = J - 1) • P(rank p (A^_i) = J - 1). 

To estimate V(iank p (AK) = J | r&nk p (AK-i) = J — 1), as in previous 
reasoning, we only have to ensure that n — J + 1 entries of the last column 
are not equal to a certain vector. On the contrary, for V{xwak p {AK) = 
J | rank p ( J 4^_i) = J) we notice that we can leave the choice of J entries 
corresponding to a non-zero minor free, but the remaining n — J entries have 
to be determined modulo p. By induction, we have 

J-2 , 1 x K-J-l 

i=o 

+ /3 n_J 0- - a(n ~ i} ) • (3 {n - J) ^ K - J - 1) ■ ( — - ) {1 + 13... p k - j - 1 ) 

= JJ(1 - • /?(»-')(*-') f— ^ J_1 (1 + (3 . . . f3 K - J ) 

i=o ^ 
which finishes the proof. □ 
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Let us consider the example of n x 2 {0, 1} matrices. We will consider 
a = (3 = \. We can construct 2 2n different matrices, 3 • (2 ra — 1) of which 
fulfill the condition that the rank is equal to 1. The probability of choosing 
at random a matrix of rank 1 is thus ^ 2 2 -i n ^ ■ The bound given by Eq. (6) 

is (1 — (^) n ) • (\) n 1 • (1 + §) which gives exactly the same value. 

The following lemma gives analogues to lemmas 5.10, 5.11 in [28] in the 
case of the ring Zy . It proves that the vectors of elements from S can also 
be treated as almost-uniformly distributed. 

Lemma A. 2. (i) Let t be a non-zero mod p vector of size n, d G 2y . 
Then the probability that a random vector x G S n is chosen such that 
t • x = d mod p l is 



V ( x : t • x = d mod p 1 ^) < — [—7] 

v J b v 



p. 

(ii) Let A G Z mxn , a matrix of rank r such that the local Smith form of 
A at p is trivial, b £ be given. Then the probability that a random 



vector x £ S r 



is chosen such that Ax 
V [x : Ax = b mod p l 



b mod p l is 



< 



s V ) 



' Ir 


" 




" R' ' 










R" 



L 



therefore transform 



Proof. For (i) the proof of 5.10 from [28] carry on. For (ii) we slightly 
modify the proof of 5.11 from [28]. Since A has a trival Smith form modulo 
p there exist two matrices L, R, det(L), det(i?) / mod p, such that A = 

where R' = [Rij] mod p is a r x m matrix. We may 

V(x :Ax = b)= V(x : R'x = [L _1 6]i.. r ) 

Since the determinant of R is non-zero modulo p there exist a r x r minor 
R\ which is non-zero modulo p. This means that we can find elements 
r u 1 ■ ■ ■ r u r °f Rii where i^ are pairwise distinct, such that r^; lk are non-zero 



modulo p. Let d 
V(x : R'x = [d}i... 



L 1 b. The probability can be further rewritten: 

= £-£ 





1 RuX! + • 


~\~ RliiX^ 


+ • 


• Rln%n 


= 31 \ 




( x h 


= (di 


V 


R21X1 + • 


■ + R2i 2 x i 2 


+ . 


■ R2n%n 


= 32 


V 


X{ 2 


= (d 2 




\ RriXi + • 


+ Rri r %i r 


+ . 


. R rn x n 


= 3r J 






= (d r 



l 2i 2 



< 



(7) 
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We use Rki k Xi k to denote that the element with index k is omitted in the 
sum. □ 

The following lemmas show that matrices of elements of S can be treated 
as almost uniformly distributed. 

Lemma A. 3. Let L,R G Z nxn be matrices such that | det (L) \ = | det (R) \ = 
1. Let I, J be any disjoint subsets of {1 . . . n} 2 (sets of index pairs). Let 
dij,(i,j) G I (resp. D s t,(s,t) G J) be any values (resp. subsets) from 
Z p i . We consider the probability of choosing a random matrix X such that 
(LXR) { j = dij for (i,j) G I under the condition that (LXR) st G D st for 
0, t) G J. We have 



V ((LXR) t] = | (LXR) st e D st ) < (ir^l 



Proof. 



V ((LXR) = d l3 A (LXR) st G D st ) 
V ((LXR h = dl] | (LXR) st G D st ) = -V v{[LXR)steDst) 

(8) 

Let V denote a set of all possible matrices [a^] such that G if 
(I, k) G J and a\\. G S otherwise. Let T> denote a set of matrices from V for 
which additionally aik = d\^ if (l,k) G X. Then Eq. (8) can be transformed 
to 



V (X G L^VR- 1 ) 

V [X G L^V'R- 1 ) 



) 

Notice, that sets P and L~ 1 T>R~ 1 (resp. D' and have the same 

number of elements. To compute the probability it suffices to count the 
number of elements in T> and T>' . The proportion is determined by the 

choice of elements from 1 and is therefore less than or equal to ^ T^l) 

□ 

The methods used to prove Lemmas A. 2 and A. 3 can applied to prove 
the following lemma. 

Lemma A. 4. Let A G Z mxn be a matrix such that the Smith form of A is 
trivial and rank(yl) = m < n. Let X, S be any disjoint subsets of {1 . . . m} 2 
(sets of index pairs). Let 6^, (i,j) G I (B st , (s,t) G S) be any values (resp. 
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subsets) from Z p i. We consider the probability of choosing a random matrix 
X such that such that [AX) i - = bij for (i,j) G X under the condition that 
{LXR) st G B st for (s,t) G S. We have 

V {(AX) l3 = b l3 | (AX) st G B st ) < (I . (9) 

Proof. Let matrices L,R = R' be as in the proof of A. 2. As in the proof of 
A. 3, we construct the sets of matrices V, V'. We have 

v((AX) -b I (AX) r P \- V ( RX& L '^) - Ih^rnRXi | L^Vj 
V [(AX) tJ - b l3 | (AX) st G B st ) - v{RXeL - lvl) ~ ^ mRX . eL -i V ,. 

where Xi denote the ith column of X and T>i(V'i), the set of all possible ith 
columns for matrices from V(V). Since (7) holds for every vector L~ 1 d of 
L^Viiresp.L^V'i), again, we can link the the probability to the number 
of elements in 1 and conclude that (9) holds. □ 

We conclude with the following lemma. 

Lemma A. 5. Let V be an k x n matrix, k < n, such that the Smith form 
of V is trival and V has a full rank. Let M be an n x k matrix with en- 
tries chosen randomly and uniformly from set S, the probability that p l < S 
divides the determinant det(VM) is at most ^. 

Proof. To check whether ord p (det (M)) > / we will consider a process of di- 
agonalization for M(0) = VM mod p l as described in Algorithm LRE of [7]. 
It consists of diagonalization and reduction steps. At the r-th diagonaliza- 
tion step, if an invertible entry is found, it is placed in the (r, r) pivot position 
and the rth column is zeroed. If no invertible entry is found, we proceed 
with a reduction step i.e. we consider the remaining (n — r + 1, n — r + 1) 
submatrix divided by p. The problem now reduces to determining whether 
ordp of an (n — r + 1, n — r + 1) matrix is greater than I — n + r — 1. 

We can consider matrix M(0) = Mq + pM\ + p 2 M2 • • • + p'M;_i, where 
matrix M k G Z" xri , k = . . . I - 2 and G Z nxn . The probability that 

an entry of M k is equal to a certain d modulo p is less that or equal T^fl > 
where is equal to \-^] by Lemma A. 4. 

In the process of diagonalization we can find matrices Lq, Rq, det (Lq) = 

det(-Ro) = 1 such that L q M Rq = diagf 1 ... 1, ... ] and L MR Q = 
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diag^.1^^1,, pLqMiRq + . . .J . Then after the reduction step we set Mq (1) = 

[(L Mii? )ij]i=r+i...n,j=r+i...n and M k (1) equal to [(LoM k+1 R )ij}i =r+ i... n j =r+ i. 
M (1) = Mq (1) + pM\ (1) + . . . and we repeat the diagonalization phase. 
By construction, the choice of Lq, . . . L k _i, Ro, ■ ■ ■ Rk-i means that certain 
entries of M are fixed and places us in the situation of Lemmas A.3,A.4. 
Thanks to that we can consider the distribution of entries of M(k) as non- 
uniform i.e. V (M(k)ij = dij mod p a \ L . . . L k -i, Ro ■ ■ ■ Rk-i) < j^TpM- 
Another way to see this is to think of the diagonalization as the mod- 
ification to <i22 in the form of 022 — I77 a i2 with an and a\i fixed by the 
previous step. Then one has one degree of freedom, say for 021 and then 022 
has to be fixed. 

We need only to consider I — 2k reductions steps as each reduction is 
performed on a matrix of order at least 2 and divides the determinant by at 
least p 2 . Since k is less than [7/2] — 1 and I < log p (S), we have N k > -p? > 

V~S and since we only consider p a < N k we have 



N k 1 p a 1 ~ p a N k ~ p a (p a + 1) p a + 1 

throughout the process. We therefore now set (5 a = ^q-j- and use it as a 
bound for (3 a (k), k = 1, 2 . . . in our calculations. 

The proof is inductive on n, the dimension of the matrix M (k) and 
I, the current exponent. We fix the diagonalization/reduction matrices 
Lq . . . L k _i, Rq . . . L k _i and consider the conditional probability V k ~ 1 = 
V (■ I Lq . . . Lfe_i, Rq . ■ ■ Rk-i)- 

First, for / = 1, [28, Thm 5.13] gives 

n 

^- 1 (pfdet(M(A ; )))<n(l-^)- 

i=i 

This transforms to 

01 



p*" 1 (p I det(M (*)))< £#<t« 



(10) 



Thus, the probability can be bounded by min ^1, and therefore by |. 

For n > 1 we will sum over all possible choices of L k and R k . We 
will divide the sum on the cases when applying L k and R k leads to the 
diagonalization of at least r entries. We call such an event E r . 
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Then forn = 2,I = 2: 
V k - 1 (p 2 | det (M (A;))) < £ P^ 1 (p 2 | (L fc M (fc) i? fe ) 22 | L fc , i? fc ) T^ 1 (L fc , fl fc ) 

4 



+ r k ^ ( P \M (k )l3 ,i, j = i, 2) < p 2 + # < ^ + (^) < 



Now we suppose inductively that (p* | det (M (k))) < 4 for all 

i < I. Then for n = 2, 2 < I < n the induction gives 



V k - X (V | det(M(fc))) < Yl ^ {p 1 \ (L k M{k)R k ) 22 \ L k ,R k )j V k ~ l (L k ,R k ) 
+ V k ~ x (p\M (k) tj ,i,j = 1, 2) P fe - X (p'- 2 | det (M (A; + 1))) < Pi (*0 + A (^) 4 ^ • 



The latter is less than 4 when 

pi 



Pi(k) 4 3p 2 <l. (11) 

With fa (k) < ^ this means that {l+l/ ^ (p+lf = (p+2 + 1/p) * < 1 which is 
fulfilled for p > 3. For primes p = 2, 3 we have to use a sharper bound for 
Pi (k). Since p' < N k and Z > 2 we have 

p ; + l+p- 1 p'- 1 + 1 p+1 

A(fc) ^ ( y + i)p VTi~ p^TT (12) 

This allows us to prove the inequality (11) for p = 3 since (I) 4 27 < 0.7. 

For p = 2 and I > 3 also (y^) 4 12 < 0.95 . For the remaining case p = 2, 
Z = 3 we can bound V k ~ x (p | det (M (fe + 1))) by 1 instead of § and then 

one can prove that p 3 (k) + pi (k) 4 < § + (|) 4 < 0.32 < |. 

Now we will consider n > 2. Again we can sum over all possible diago- 
nalization and reduction steps combinations and the resulting bound for the 
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probability is 

V k - 1 (> | det(M(fc))) < V k ' x (p\M (k) tj y h3 < n ) 



n—l 



+ E E ^ (P\ ( L « M R k)ijVi j<n-r I £fc, flfc) (£fc, Afc) 

r=l L k ,R k eE r 
n-2 

+ E E pfe_1 (pi ( LfcM v ^<™- 1 L *> • 

r=n— l+l L k ,R k £E r 

V k - 1 (L k , R k ) V k (p l ~ n+r \ det (M(k + 1))) 

n 

L k ,R k &E n -i i=l 
n-2 

+ E E /3i(fc) ( ^ r) V fc - 1 (L fc , J R fe )P fe (^-"+ndet(M(fc + l)))+A(fc) 

r=n— (+1 L k ,R k £E r 

(13) 

for Z < n and similarly for Z > n 

P^ 1 (V | det (M (&))) < P^ 1 (p|M {k) tj Vij< n ) P^ 1 (p'- n | det (M(fc + 1))) 

n-2 

+ E E (pi v ^<«- 1 L ^ ^) • 

r=l L k ,R k £E r 

V k - 1 (L fc) P fc det (M(fc + 1))) 

+ E ^ (P l \ ( LkM ( fc ) I ^ ( L *> ^) 

L k ,R k (LE 

n — l 

< A (fc)" 2 V k ~ x (p l ~ n \ det (M (fc + 1))) + A (*) 

n-2 

+ E E a ( fc ) (n_r)2 ^ ( L *> R ") v k ~ l {p l ~ n+r \ det ( M ( k + !)) i fl*) • 

r=l L k ,R k £E r 

(14) 

Again, we can use the induction to get the bound V k (p 1 " 1 | det (M(fe + 1))) < 
Then, we can then bound both sums by 

V k - X (p l | det (M (k)j) < £ ft (^r 2 4?+A (fc) < ^#^7 ^^T+A (*o 

(15) 

38 



To prove the inequality V^ 1 (p l | det (M (k))) < 4, we have to consider 
several cases. For p > 3 we use the bound Pi and Pi for Pi (k) and $ (A;) 
respectively. Then we have 

3-2 4 p 2 (p + l) 2 2 1 48p 2 (p+l) 2 1 48p 2 



p»((p+l) 5 -p25) ^ + 1 P l P l (p+lf-( P +l)2* P l P l {p+lf 

2 1 48p 2 2 1 48p 2 3 

< ^7 + y 25^2 + 4 . 5p 2 + 6 . p 2 + 4 . 5 + ! _ 16 < ~i + Ji ^2 < ^7" 

For p = 3 it can be explicitly checked that V k ~ x (p l | det (M)) < ^ us- 
ing the bound for p\ (/c) (notice that N k > p'). In this case we get 

1 3 (f) 432 I 2 ,lo 7 r 

3^ (l-(3|) 5 ) + 37 < ^ 

For p = 2 we have to consider 2 2 , 2 3 , 2 4 and 2 l for Z > 4 separately and 
use the sharper bound from Eq. (12). Let us rewrite (13) and (14) in this 
cases. 

• 1 = 2: 

V"- 1 (2 2 | det (M (k))) < Pi (kf+(3 2 (k) < Pi {kf ] +p 2 (k) . 

7^2 1 - ft ( fc ) 

As A (A;) < |Jt we have 0.65 < 0.75. 

• / = 3: 

n 

V k ~ X (2 3 | det (M (fc))) < £ Pi (kf + Pi {kf ■ 1 + p 3 (k) 

< Pi {kf 1 + /?i(A;) 4 + 
1 - Pi (k) 



2 4 



As p x (k) < §±± we have 0.33 < 0.375. 
• / = 4: 



V k ~ x (2 4 | det (M(fc))) < Pi(k) l2 + Pi{kf- 1 + P\{kfV k (2 2 | det (M{k + 1))) + 



i=4 



" /3l(fc)16 l-ft(A:)» + A(A;)9 + A(fc)4 I + /?4(fc) - 

As p x (k) < ^ we have 0.18 < 0.1875. 
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• / > 4: 

We use inequality (15) with ft (k) bounded by *£±i. We get T^" 1 (2 ? | det (M (A;))) 
is less than £ ( ^^ff^l^^ + 2) < 2.92^ < § . 

We have thus proven that V k ~ l (p l \ det(M(fc))) < ^ for every / > and 
every size n of M(k). Thus, V(p l \ det(VM)) is also less than or equal 4- □ 
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